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The Quantum Monte Carlo simulation of the two-dimensional Emery model of Cu02 plane of 
hight Tc oxide superconductors were performed. The method based on the direct-space proposed 
by Suzuki and al., Hirsch and al. was used. Contrary to the method based on the Hubbard- 
Stratonovich transformation, the states generated by this method are basis states in occupation 
number representation, i. e. configurations of fermions can be observed on real two-dimensional 
array. Energy and specific heat were computed for different dopings during decreasing temperature. 
The specific heat curves show peaks at low temperature which could be assigned to electronic 
transitions. Quantity similar to current-current correlation function were computed. The static 
electric conductivity curves obtained by this way show metal-insulator transitions for doping S — 
and doping 5 = 1, and two different metallic behaviours for intermediate dopings. On the direct- 
space states generated at low temperature and doping 5 = 0, the fermions form antiferromagnetic 
loops while they form antiferromagnetic chains for other dopings. The loops can be related to 
circulating currents and the chains to the stripes. The antiferromagnetic loops seem to appear 
when the conductivity becomes zero while the conductivity increase with the numbers of chains but 
superconductivity is not unambiguously evident. 

PACS numbers: 71.27.-|-a, 71.10.Fd, 71.30.-fh, 74.20.Mn 



I. INTRODUCTION 

The Emery model[l| of the Cu02 plane was proposed 
and extensively studied to explain the behaviour of hight 
Tc oxide superconductors. This model seems to repro- 
duce the essential phenonema observed on materials al- 
though the evidence of superconductivity is not totally 
proven for the repulsive Emery model. Many results are 
obtained by numerical simulations using different quan- 
tum Monte Carlo methods. These methods are based 
on the Hubbard-Stratonovich transformation 2] 3] which 
does not generate actual fermions configurations. 
In this paper, we present results on two-dimensional re- 
pulsive Emery model obtained with a method based on 
the direct-space proposed by Suzuki and al. [1| ^ and 
Hirsch and al.QQ- At fixed temperature, this method 
allows to generate some of the most representative occu- 
pation number basis states of the model. These states are 
used to compute average values of energy, specific heat 
and rough static electric conductivity. 
The paper is organized as follows: 

• in section 2 the Emery model is recalled, 

• in section 3 the numerical method is described, 

• in section 4 the validity of the method is tested with 
the two-dimensional Hubbard model, 

• in section 5 the energy and specific heat curves are 
shown, and some direct-space states are presented, 

• in section 6 we define a way for computing a rought 
value of conductivity and present the computed 
curves. 



• in section 7 we analyze the results, 

• in section 8 some concluding remarks are given. 

II. EMERY MODEL 

The studied system is the Cu02 plane which have the 
structure shown in Fig. [T] One considers the behaviour of 
holes on this array and uses the occupation number rep- 
resentation. The single-particle states used to construct 
the basis states of the representation are the hole Wannier 
states localised on each site. According to atomic orbitals 
involved for copper and oxygen atoms, the Wannier state 
on all copper sites are labeled d while Wannier states on 
all oxygen sites are labeled p. It will be assumed that 
the Hamiltonian of the system is an extended Hubbard 
Hamiltonian 

H ^ tdp {d]^^Pi,a + hcj +Ud'Y nd,iind,i'[ 
+Up ^ rip^iiUp^i-i + Ed^^ nd,i + Ep^^ rip^i 

i i i 

nd.iripj (1) 

The indice (i, a) labels the hole Wannier state localised 
on the site i with spin a. The operator d| ^ creates a hole 
in the copper Wannier state {i, a) while the operator di u 
destroys a hole in this state (i, tr). The operator p| ^ cre- 
ates a hole in the oxygen Wannier state (i, a) while the 
operator pi^^ destroys a hole in this state {i, cr).The oper- 
ator nd,i,a = d] „di^a is the occupation number operator 
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FIG. 1: The structure of the Cu02 plane. The squares indi- 
cate Cu sites and the stars oxygen sites. 



the presence of the off-diagonal hopping interaction, the 
eigenstates of the Hamiltonian, Eq.(IT|), are not basis 
states, in such a way that the computation of the traces 
in Eqs.(I2) and ^ is problematic. The method proposed 
by Suzuki and Hirsch allows to get round these difficul- 
ties 130. 

The Hamiltonian of the system is decomposed in sev- 
eral sub-hamiltonians Hr comprising the creation and 
annihilation operators of only certain sites. Because of 
the fermion commutation relations some of these sub- 
hamiltonians do not commute. This breakup is not to- 
tally arbitrary, the sub-hamiltonians are chosen like each 
sub-hamiltonian cans be decomposed again in several 
non-overlapping sub-system Hamiltonians K^^k commut- 
ing. 



of the copper state {i,a) and Ud^i = nd,ii + is the 
holes number operator on the copper site label i. np,i^cr 
and np,i are the same operators for the oxygen states. 

indicate that the summation is performed on the 
nearest neighboring sites. 

The values of the interaction parameters tdp^ Ud, Up, Ed, 
Ep and Vdp were estimated by differents authors 3 [9|] ilO|] 
but the results of our simulation being very sensitive 
these values we choose them after several tests. The val- 
ues are chosen such that the low-temperature conductiv- 
ity variations are large enough. The retained values of 
the interaction parameters are, in eV : the hopping ma- 
trix element tdp — —1.5, the onsite repulsive Coulomb 
energies C/^ = 9 and Up ~ 0, the site energies Ep — i 
and Ed = 0, and the Cu-0 intersite Coulomb energy 
Vdp = 0.75 Figdl 

In this model, the vacuum corresponds to the state in 
which all the oxygen p orbitals and copper d orbitals are 
fully occupied. There is one hole by elementary cell for 
the undoped plane. 



III. NUMERICAL METHOD 



A. Simulation method principle 



H ^ [Hr,Hr']^0 T ^ r' (5) 

T = \ 

raj. 



k=l 



The previous breakups are realized with the aim of using 
the Trotter formula, Eq.([7]), which allows to get round 
the non-commutativity problem of the sub-hamiltonians 



lim 

n — *oo 



P r 



n 



exp Hr 

n 



(7) 



Using the Trotter formula the partition function Z reads 
Z = lim Zn (8) 

n — >oo 

where Z„ is a partition function approximant 



Zn = tr ■ 



n 



exp Hr 

n 



(9) 



Inserting np complete set of states between operators the 
approximant Z„ is given by 



The principle of the simulation method is detailed in 
references Q [H| • We recall here the essential caracter- 
istics of the method. Within the canonical ensemble the 
average value of an observable, O, is given by 



(O) = tr (DO) 



(2) 



where D is the density operator and Z the partition func- 
tion 



-I3H 



D 



Z = tr {e-P") 



Z 

I3H\ 



(3) 
(4) 



Here, H is the Hamiltonian of the system and (3 = l/ksT 
the inverse temperature. In the present case, due to 



Zr^ = (*0|exp (--Hp) \^np-l) 

(5'„p_i|exp [~-Hp^i ) |1'„p_2 



(^-ilexpl -^i/i ) l^-o) (10) 



Here represents the configuration of the np states 
and can be seen as the state of a (d -1- 1) dimensional 
classical system, d being the dimension of the studied 
real quantum system. indicates that the sum runs 

over all the possible configurations. This is equivalent 
to divide the imaginary-time interval ^ r ^ /? into n 
slices of width Ar = (3/n. 
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An approximant of the energy of the system is 
Substituting Z„ from Eg. pop . we obtain: 



(11) 



Un = Pn{['^^])En{[^^]) (12) 

{[*.]} 

En ([*.]) = E J (13) 

,=0 (vf.+il exp 

^«([*z]) = ^(*olexpf-^i/j,) |*„p_i} 
Z„ \ n / 

...(I'llexp (^-^i/i^ I'J'o) (14) 

The sub-hamiltonian index r of each matrix element is 
a function of the state index j, so r = 1 + {j modp), 
= |*o) and 



E ^"([*^]) = i 

{[*.]} 



(15) 



We can consider that each configuration [^i] of the 
{d+ l)system has an energy i?„ {[^i]) and a probabil- 
ity factor Pn {[^i])- The computation of the average en- 
ergy of the system is realized with a usual Monte Carlo 
method like Metropolis algorithm. 
Using Eq.([6]), the computed value of Un is: 



1 



(16) 



i=l r=l fe=l 



(vl/,+i|exp(~^E^r,i)|*, 



(17) 



where Np is the number of kept configurations, a is the 
index of configurations replacing [^i] and i is the slice 
index. The indices i,j,r are linked by the numbering of 
the np states and verify j = {i — l)p + r — 1. 
The breakup of sub-hamiltonians Hr allows an important 
simplification which necessitates another approximation. 
All the sub-system hamiltonians Kj..i of the same sub- 
hamiltonian Hj. commute and each hamiltonian K^ i acts 
only on the state of one sub-system, thus, one writes 



I exp 



f3 



E 



1=1 



j+i,i\exp 



13 



-Kr,l ) \ifj,l) 

n ' 



(18) 



where the state \<fj^i) is the state of the I sub-system of 
the r sub-hamiltonian of the i slice. The symbol < — > 



means that the left expression is replaced by the right 
one. This operation means that the space of the states 
of the system is considered the tensorial product of the 
spaces of the sub-systems states. It is not true for a 
fermion system, because of the state antisymmetry, the 
creation and annihilation operators are defined in the 
state space of the whole system. There is not strictly 
equality. Implicitly, this means that occupation numbers 
of the others single-particle states are not taken into ac- 
count. This approximation is used by all the authors who 
present results obtained with world lines simulations. 
The same factorization for the numerator of Eq. p7|) is 
applied, finally, the computed values are: 



U'n - {E'n^a) 



(19) 
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1 {ipj+i^k\Kr,keyiY>{-^Kr,kj Wjm) 



n 



EEE 



r=l k=l (-^j-l-l.fcl exp (-fi^r.fc) W].k) 



(20) 



1 f (3 \ 

Pn,a = ^ n n n (V'i+l.'cl exp l --Kr.k j 
" i=l r=l k=l ^ ^ 



(21) 



n p m-r 



^« = E n n n (^i+i.^! --Kr,k i \^j,k) 



i—l r—1 k—1 



(22) 



Some probability factors P„ {[^i]) or are negative: 
it is the "sign problem^^ . Hirsch and al. got round this 
difficulty by computing (_E„sgnP„) @. 
The specific heat is calculated with the formula 



dU 
'dp 



(23) 



where k is the Boltzman constant. 

Using a similar method that this used for we obtain 
the fluctuation formula 



C 
R 



p'{{E'n,.-Ky) 



(24) 



where C is the molar specific heat and R the perfect gas 
constant. 



B. Simulation parameters 

Our simulation model of the Cu02 plane contains 6x6 
elementary cells with periodic boundary conditions. Fig. 
[2] shows the breakup of the model. There is one sub- 
system around each copper site. All sub-systems are 
identical and are composed by one copper site and four 
oxygen sites. These sub-systems are grouped together in 
two sub-hamiltonians {p = 2) in such a way that the sub- 
systems of one sub-hamiltonian do not have common site. 
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The sub-system hamiltonians are identical and called K. 
Using the numbering 1 of the sub-system sites given in 
the figure [3l the sub-system hamiltonian reads 

K= tdp (4cr {Pla + P2cr + Via + P5a) + hc) 

a 

+Ud ridsindsi: 
Up , 

('^pU'^piT + np2inp2^ + np4inp4^ + np^iUp^^) 
Ep 

+Ed nd3 + [ripi + np2 + Upi -|- rips) 

+Vdp rids {npi + np2 + Up^ + Up^) (25) 

The factor 1/2 is introduced to take into account that 
each oxygen site, numbered 1,2,4 and 5, belongs to two 
sub-systems. 

The state space dimension of each sub-system being 2^" 
the method necessitates the diagonalization of only one 
(1024 X 1024) matrix. A constant value is added to the 
eigenvalues in such a way that the ground level is zero. 
The imaginary-time interval is divided into twenty four 
slices (n — 24). Thus the (2 -1-1) dimensional system is 
composed of forty eight 2D-lattices [np = 48). 
Each simulation begins with thermalization. Then 
decreasing temperature (100 points) is programmed, at 
each temperature point the averages are computed on 
100 configurations of the 3D-array. There is about 50,000 
Monte Carlo steps between two kept configurations. 
Simulations were performed for different dopings 5. For 
6x6 elementary cells, there are 36 holes (18 ] +18 |) 
for doping (5 = and 72 holes (36 T +36 j) for (5 = 1. 
The energy average U = {Ei) of the model is calculated 
and the ratio C / R is computed with the energy fluctua- 
tions formula, Eq. (I24p . 

The comparison of C/R with the derivative of the cell 
energy can be used as criterion of the numerical method 
quality. 




FIG. 2: (Color online). The breakup of the Cu02 plane into 
two sub-hamiltonians with periodic boundary condition. The 
grey (or red) sub-systems belong to the sub-hamiltonian 1, 
and the dark sub-systems belong to the sub-hamiltonian 2. 

off-diagonal matrix elements {Lpj+i\K\ipj) is calculated 
for three different numberings of the sites. The figure [3] 
shows the corresponding situations. The element values 
are 

numbering 1 

{^,+4\K\^j) - (0, 1,1,0, m\h 1, 0, 0, 0) = -tdp ^^^^ 
numbering 2 

(27) 

(<^j+i|if|</jj) = (0,1, 1,0, 0|7^| 1,0, 1,0,0) = trfp ^ ' 



C. Influence of the state numbering 

Considering that the sub-system hamiltonian K does 
not modify the particle numbers, the matrix of the hamil- 
tonian K cans be decomposed in sub-matrices corre- 
sponding to sub-spaces of fixed particle numbers. The 
maximum dimension of these sub-spaces is 100. It corre- 
sponds to the four sub-spaces of 2 spins up and 2 spins 
down, or 2 spins up and 3 spins, or 3 spins up and 2 
spins down, or 3 spins up and 3 spins down. The dimen- 
sion of the sub-space corresponding to 2 spins up and 
zero spin down is C| = 10. This sub-space is chosen to 
test the influence of the single-particle state numbering 
on the matrix of the hamiltonian K and consequently on 
the probability factor sign. 

The numbering principle is chosen such that the numbers 
of the single-particle states correspond to the sites num- 
bers in the present case where there is only spins up. An 



numbering 3 

{ip.j+i\K\^j) = (l,l,0,0,0|i^|0,l,0,0,l) = -tdp ^^^^ 

The sign of the off-diagonal matrix elements are depen- 
dent of the state numbering. Indeed the change-of-basis 
matrices are unitary (orthogonal) matrices and the ele- 
ments are or ±1 with only one +1 or one —1 by row 
and column, the others elements being 0. They are per- 
mutation matrices except that some matrix elements are 
— 1 instead +1 because of the antisymmetry. 
These remarks presented for the particular (10 x 10) sub- 
matrix hold for the all sub-matrices of the hamiltonian 
K. For a numbering change, the change-of-basis matrices 
(and the inverse matrices) act on the sub-matrices of the 
operators exp(— /^if) and if exp (— ). Thus the sign 
of some elements of these sub-matrices change when as 
the absolute values are equal. For a fixed particle num- 
ber sub-space, the modified elements of the two matrices 
are the same. Thus the sign of each ratio of the sum in 
the relation (|^D|) is not change, whereas the sign of some 
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D. The probability factor signs can be ignored 



FIG. 3: Three difFerent numberings of the sub-system sites. 
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FIG. 4: (Color online). Energies of 3D array configurations 
versus temperature. The probability factor signs correspond 
respectively to the three numberings of the figure [3] 



factors of the product in the relation ([2T|) change. 
These remarks show that the energy E'^ ^ , Eq. ([20| , of 
each configuration of the 3D array is independent of the 
numbering while the sign of the probability factor ^ , 
Eq. pT|) . cans change. The figure |4] shows the variations 
of the energies E'^ ^ of different 3D array configurations 
with the temperature, and the probability factor signs, 
for the three numberings of the figure [31 



The partition fonction Eq. ((22)) . and the proba- 
bility factor Pj' ^ , Eq. (I^T|) , can be written 



ZL 



pi _ 



n p m-r 



(29) 
(30) 



i—l r—1 k—1 



(3 



n ' 



(31) 



where Oq, is relative to 3D array configuration indiced a. 
The energy E'^^^ and the absolute value of Oq, arc prop- 
erties of the 3D array configuration, they depend only on 
the configuration whereas the sign of Oq, is dependent of 
the state numbering. 

During numbering change, some positive terms aa be- 
come negative and some negative terms become positive. 
The predictions of the method must be independent of 
the state numbering, the partition funcion Z'^ and the 
energy average must be unchanged, thus 



O-aEa — — ^ aaEa — 



E 



(32) 
aa\Ea (33) 



where Ea replaces -E'j {an-,p} represents the set of 3D 
array configurations of which the terms become pos- 
itive and {ap^n} the set of 3D array configurations of 
which the terms Oq become negative. The energy aver- 
ages of these sets of configurations are equal: 



1 a^Ea 



}' 



(34) 



}' 



This result is true for all the different possible number- 
ings of the single-particle states of the sub-system. 
The factor Ua of each configuration a is positive or nega- 
tive depending on the numbering and the energy of each 
3D array configuration is a contributation to averages 
computed with positive and negative factors. Thus we 
can deduce that, for a given numbering, the energy av- 
erages computed separately for the positive en negative 
probability factors are equal: 



{a<o} 



(35) 



where {a>o} and {a<o} represent respectively the sets of 
3D array configurations with positive and negative prob- 
ability factors. This equality between the energy averages 



computed separately with the negative'^ and ''''positive'^ 
configurations was aheady remarked Ours compu- 
tational results confirm this equality. This behaviour is 
independent of the number of kept configurations, the 
number of slices, the size of the model and the tempera- 
ture. 

Using a well-known arithmetic result one cans write 



E{a} \0-a\Ea 



I]{q} 



(36) 



S{a>o}""+E{a<o}«o 



(37) 



where Pq, replaces ^ . The equality of terms and 
([57)1 indicates that the sign of the probability factor cans 
be ignored. We use this remark for our simulations. 
The computational results are actually independent of 
the single-particle state numbering. 



IV. TEST OF THE VALIDITY OF THE 
METHOD 

The validity of the method is tested with the two- 
dimensional Hubbard model on a square lattice at half 
filling. This model was extensively studied with different 
methods [i2|[3[3- The model and simulation param- 
eters were chosen to compare our results with the results 
of ref. (fig. 6 and fig. 8). The Hubbard hamiltonian 



H = t 



E (4. 



he 



UY^n 



(38) 



The imaginary-time interval is divided into eight slices 
(n = 8), the lattice is 6 x 6 with periodic boundary con- 
ditions, with eighteen spins up and eighteen spins down 
(18 t +18 i). The interaction parameters are t = —1 
and U — 2 or U — 10. All sub-systems are identical 
and are composed by four sites. These sub-systems are 
grouped together in two sub-hamiltonians {p — 2). 

The Fig. [5] displays the molar specific heat curve as 
a function of temperature kT for U = 2. The curve 
is an average of eighteen simulations while the curve 
for [/ = is an average of six simulations. There are 
100 temperature points by curve (20 points by decade). 
For the sake of clarity all the error bars are not shown. 
These curves are very similar to the curves in ref. 



t=-l; U=2; 18?+ 18J. 

- - t=-l; U=0; 18?+ 18J. 
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FIG. 5: (Color online) Molar specific heat versus temperature 
for the Hubbard model at half filling, t — —1, U — 2 and the 
lattice is 6 X 6. The dashed line curves (blue) correspond to 
17 = 0. 



[131 . In particular, a two peaks structure is observed, 
with a broad high temperature peak and a narrower 
peak at low temperature. However, the temperature 
range is translated by a factor four toward the low 
temperature and there is a ratio four between the C/R 
amplitudes. This ratio is consistent with the definition 
C = dU/dT. Indeed, in ref. [l2| the slice number n 
is not constant while n = 8 in our simulations, this 
difference can explain the temperature translation. The 
intensity of the low temperature peak of our simulation 
is larger than in ref. [l^l- Note that the error bars are 
large in the temperature range of the low temperature 
peak. This indicates that there are large fluctuations 
in this temperature range (0.033-0.05). A simulation 
with 200 temperature points in the range (0.01-0.1) was 
performed to understand this problem. Fig [5] shows the 
energy curves in this temperature range. The dashed 
line curve with error bar corresponds to the average of 
eighteen simulations with 100 points in the temperature 
range (0.001-100). The solid line curve is the result 
of the simulation with 200 points in the temperature 
range (0.01-0.1). The energy of the states generated into 
the temperature range (0.033-0.05), oscillate between 
two values. The error bars are very small outside this 
temperature range. This indicates that a first order 
transition occurs. This explains the large fluctuations 
of the low temperature specific heat peak. The results 
in ref. [lj| are obtained with the determinantal QMC 
simulations, based also on the Trotter formula. The spe- 
cific heat curves C (T) are evaluated by differentiation 
of the energy data with about 10 temperature points 
by decade. Our specific heat results are raw data of the 
simulations, computed with fluctuation formula. These 
differences of methods and numbers of points explain 
the differences between the results. 

The Fig. [7] displays the molar specific heat curve for 
U = \Q (average over 18 simulations). This curve is very 




Temperature 

FIG. 6: (Color online) Site energy versus temperature for the 
Hubbard model at half filling, t — —1, U — 2 and the lattice 
is 6 X 6. The dashed line curve (black) corresponds to the 
average of 18 simulations with 100 points. The error bars 
are very small outside the transition range. The solid line 
curve (red) is obtained with 200 points in temperature range 
(0.01-0.1). The dot line curves (blue) are extrapolations of 
the energy curves. 



similar to the curve in ref . [12] ■ The intensity of the peak 
for t = is slightly higher than the high temperature 
peak but the uncertainty can explain this difference. 
For t — the sub-system hamiltonian is diagonal and 
there is no sign problem. The error bars indicate that, 
for t = —1 and U = 10, the algorithm generates states 
which stay a long time in local minima for increasing and 
decreasing temperatures. This problem is less important 
for the Cu02 plane. Thus a small number of simulations 
is required to obtained small error bars. 

The test of the validity of the method is completed with 
simulations for [/ = (fig®. For [/ = 0, the half fill- 
ing (18 t +18 i) corresponds to two no-correlated quan- 
tum gases with the same energy and specific heat. The 
simulations for 18 t and no spin down give identical re- 
sults for the spin energy and the spin specific heat. The 
ground state energy by spin is 2t. For only one spin (If) 
on the array the antisymmetry of states is not relevant, 
thus there is no sign problem. The spin energy and spe- 
cific heat for one spin are twice the values computed for 
18 t +18 ], or 18 t +0 i- As expected, the ground state 
energy is At. The value of the ratios of spin energy and 
specific heat computed for 1 t -1-0 | and 18 t +18 i indi- 
cate that the sign problem is correctly got round in our 
simulations. 

The similarity between the results of ref. [l^ and our 
results, and the good overall agreement between simula- 
tion results with no sign problem and simulation results 
with sign problem confirm that the sign is not an actual 
problem for our simulations. The simulation method is 
valid. 
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FIG. 7: (Color online) Molar specific heat versus temperature 
for the Hubbard model at half filling, t = —1, U = 10 and 
the lattice is 6 x 6. The dashed line curves (blue) correspond 
to t = 0. 



V. CUO2 PLANE: FIRST RESULTS 

The Figs. |9] and \W\ display the energy and specific 
heat curves versus temperature kT (in eV) for the Cu02 
plane. The curves of Fig. [10] are the averages of five 
simulations. Whatever the algorithm quality, one no- 
tices that the low temperature state of the model is the 
ground state. We observe a crossing-point on the set of 
specific heat curves with C/R ~ 2.7. This phenomenon 
is often observed for in correlated systems [I^- Fig. [11] 
shows one of the 48 2D-lattice states of the 3D-array for 
doping S — 0, while Fig. [T^ shows one 2D-lattice state 
for S = 0.33 (24 t +24 j). These 3D-arrays are the last 
generated at the lowest temperature {kT ~ 2.510"''). 
For (5 = 0, we remark some antiferromagnetic loops, and 
antiferromagnetic chains for 6 — 0.33. Only the oxygen 
sites are occupied in AF loops. It is the same for AF 
chains except that the end chain sites can be copper 
sites. The loops can be related to the circula ting current 
proposed by C. M. Varma [3 , [i3 , [ill , [3, , il 
while the chains recall the stripes proposed by several 
authors [22] 23^|. This comforts us to study the average 
temperature dependence of the chain and loop numbers. 
Fig. [13] shows these temperature dependence. It appears 
clearly that, for 6 = 0, the C/R peak is associated with 
the chain and loop numbers temperature dependence. 
There is no evidence for the other dopings. 

Fig. [13] displays one 2D-lattice at temperature kT ~ 
0.025 for S = 0. A tendency to antiferromagnetic configu- 
ration, where only copper sites are occupied, is observed. 

For (5 = 1 all the low temperature 2D-lattices are iden- 
tical to the 2D-lattice of Fig. [T5] 

The numerical method quality is tested by compar- 
ing the specific heat computed with Eq. ([24]) with the 
derivative of the energy curve for doping (5 = 0. Running 
average is applied to the energy curve before derivative. 
The derivative must be divided by the number of elemen- 




FIG. 8: (Color online) Spin energy and spin specific heat ver- 
sus temperature for the Hubbard model for different fillings. 
t = —1, U = and the lattice is 6 x 6. The solid line curves 
(blue) correspond to one spin up on the array. 

tary cells to obtain C/R. The Fig. [TBI shows the curves 
obtained with the two methods. There is a good agree- 
ment for low temperatures. 



VI. CONDUCTIVITY 

The conductivity of the model for the different doping 
is an important question. The present simulation method 
does not aUow rigorous conductivity computation be- 
cause of current-density operator defined by Scalapino 
and al. ^24*1 is not diagonal and can not break up into 
sub-system operator [7] . To obtain rough approximation 
of the current-density at imaginary-time r = ^"^j we 
calculated: 

jx {I; m) = ^ ric^l-m (1 - "<j,i';m) na,l':m+p 

(7 

(1 — n.crd-m+p) — ricf^V-m (1 — n^,l;m) 
ncr,l;m+p (1 — ?^cr,i';m+p) (39) 

The coordinates of the I and I' sites verify x// = xi + 1 and 
yi' = yi. na,i;m is the occupation number of the single- 
particule state \a,l) of the 2D-lattice numbered m. This 



FIG. 9: (Color online). Energy and specific heat versus tem- 
perature for different dopings of the Cu02 plane. 

expression reads 

{lpyn+p\ ^ C t<x,i' Ca.l - C t<x,i Ccr,i/|'0m) (40) 
a 

where the state \ipm) corresponds to the state of I and V 
sites. The computed conductivity is: 

^ np Ns 

^xx (m) = — jx {i; m + i) (I; i) (41) 

^ 1=1 1=1 

where Ng is the sites number. The DC conductivity is 
obtained after a discrete Fourier transform. Figs. [I7l 
and [m show the DC conductivity for the different dop- 
ings. If this calculation is relevant, we observe two metal- 
insulator transitions for ^ = and i5 = 1. 



VII. DISCUSSION 

A. Doping 5 = 

By contrast with the behaviour usually expected, the 
low-tcmperature state is not antiferromagnetic (AF) with 
only occupied copper sites, this can be observed on the 
figure [TT] As the temperature is decreased, the system 
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FIG. 10: (Color online). Energy and specific heat versus tem- 
perature for diflFerent dopings. 



state seems to go to AF state (Fig. [H]) and the conduc- 
tivity decreases but it rapidly increases before dropping 
to zero. The ground state is essentially composed with 
antiferromagnetic loops. Apparently this ground state is 
a new insulator phase which is not a Mott insulator as the 
AF state. The sharp conductivity increasing is difficult to 
explain, while the conductivity dropping corresponds to 
the transformation of chains into loops. Indeed the tem- 
perature dependence of the chains and loops numbers 
have similar dependencies that the conductivity. The 
specific heat peak corresponds to this transformation. 



FIG. 11: Typical 2D-lattice for 5 = at kT = 2.5 10"* eV. 
Points indicate non-occupied copper sites. Sign -I- indicates 
spin up and sign - spin down. The dot lines show the antifer- 
romagnetic loops. 
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FIG. 12: Typical 2D-lattice for 5 = 0.33 at kT = 2.5 lO"* eV . 
Points indicate non-occupied copper sites. Sign -I- indicates 
spin up and sign - spin down. The dot lines show the antifer- 
romagnetic chains. 



B. Other dopings 

1. Conductivity and chains 

The conductivity curves for 6 — 0.16, 5 ~ 0.33 
and (5 = 0.5 show that the system exhibits a metallic 
behaviour at high and low temperatures with a semi- 
conductor like behaviour for the intermediate range. 
Thus, these two metallic behaviours must be different. 
At low temperature, the conductivity shows a maximum 
for 5 = 0.33. The low temperature metallic behaviour 
seems to be associated with the chain number but the 
temperature dependence of the conductivity and chain 



number curves are not similar. A more detailed analysis 
of the chains is needed. 

There are three different types of chains: chains without 
"order" , antiferromagnetic chains and antiferromagnetic 
chains with the two end sites being copper sites. The 
last chain type is composed of three kinds of chains 
according to whether the spins of the holes of the chain 
ends are identical or opposite. The numbers of this 
different chains are computed. Fig. [12] shows the curve 
of the number of chains with opposite end spins. These 
chains, called TO, have an even holes number. Fig. [12] 
displays also the curve of number of chains with spin up 
ends, which have an odd hole number. These chains are 
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FIG. 13: (Color online). Chain and loop numbers versus 
temperature for different dopings. 
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FIG. 16: (Color online). Comparaison of C/R curves for S = 
0. Solid line curve corresponds to fluctuations computation, 
dashed line curve is obtained by energy derivative. 
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FIG. 14: Typical 2D-lattice for 5 = at kT ~ 0.025 eV. 
Points indicate non-occupied copper sites. Sign + indicates 
spin up and sign - spin down. 



is associated vifith the TO chains while the high tem- 
perature metallic behaviour is independent of the chain 
number confirming the difference between these two 
metallic behaviours. This high temperature behaviour 
probably corresponds to conventional metal, i. e. Fermi 
liquid. The comparison of TO and Tl chain numbers 
for dopings S — 0.33 and S — 0.66 shows that the Tl 
(and T-1) chains do not contribute significantly to the 
conductivity. By contrast, the Tl (and T-1) chain 
number increases with the semiconductor like behaviour 
as temperature decreases. 



called Tl. The results for the last type of chains, called 
T-1, are similar to Tl and are not shown. 

The temperature dependence of the TO chain number 
and the conductivity are similar. This relation between 
the conductivity and the TO chain number seems to 
indicate that the low temperature metallic behaviour 



2. Specific heat 

For S = 0.16, S — 0.33 and S = 0.5 the specific heat 
peaks are associated with the TO chain formation unlike 
the peak for 6 = which is correlated to loop formation. 
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FIG. 18: (Color online). Conductivity versus temperature for 
different dopings. 



However, all the peaks appear at the same temperature 
as if they originate from the same phenomenon. 
Fig. [^nishows the number of chains which are not antifer- 
romagnetic and the number of antiferromagnetic chains 
which do not have two copper end sites, i. e. which are 
not TO, Tl or T-1 chains, these chains are called Ca-t 
chains. The peaks of the Ca-t chain number derivatives 
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FIG. 19: (Color online). TO and Tl chain numbers versus 
temperature for different dopings. 



(not drawn) are centered at the same temperature that 
the specific heat peaks with a similar shape. We can de- 
duce that these specific heat peaks correspond, partly, 
to transformation of Ca-t antiferromagnetic chains into 
antiferromagnetic loops for 6 = 0, and TO antiferromag- 
netic chains for the intermediate dopings, during cooling. 
The Ca-t chains are mainly composed by antiferromag- 
netic chains with one end site being copper site. Indeed 
the ratio of the number of the chains without copper site 
end and the number of chains with one copper site end 
is smaller than 10^^. Thus, we can consider that the 
specific heat peaks correspond to "capture" of one hole 
on the second copper site end of each antiferromagnetic 
chain with one copper site end already occupied by a hole 
of opposite spin. The resulting TO chains have an even 
number of holes. 

The specific heat peaks can be cornpared to the oxide 
superconductor experimental data [25j.[26j which show 
that the height of the specific heat peaks decrease with 
the doping when 6 varies from 0.03 to 0.84. Our compu- 
tational results are in good qualitative agreement with 
these experimental data. 
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FIG. 20: (Color online). No-antiferromagnetic and Ca-t 
chain numbers versus temperature for different dopings. 



3. Superconductivity 

Given the sudden conductivity increasing, the super- 
conducting transition is questionable and the role of the 
TO chains remains an important question. The maximum 
low temperature conductivity is not upper than the max- 
imum high temperature conductivity {kT ~ 0.075 eV). 
The ratio of the maximum of conductivity at low tem- 
perature and the minimum is about three. However, the 
sharp conductivity increasing can point out a tendency 
to superconductivity which is related to the TO chains. 
The coherence of all the results supports the fact that the 
chosen approximation of conductivity A^x is relevant. 



VIII. CONCLUSION 

This study confirms that the Emery model exhibits a 
complicated behaviour at low temperature and probably 
captures the essential physics of some materials with 
Cu02 plane. However, no sure evidence of supercon- 
ductivity is found for the interaction parameters used 
in our simulations. Our results show metal-insulator 
transitions for S — and 6=1 dopings. For intermediate 
dopings two different metallic behaviours are evidenced. 
A more detailed study shows that the low temperature 



metallic behaviour is due to antiferromagnetic loops and 
chains. The physic of these antiferromagnetic chains 
(or stripes?) seems dominate at low temperature. The 
antiferromagnetic loops can be related to the circulating 
currents. 

As for all simulations in general, we must point out 
the weakness of this method. The periodic boundary 
conditions combined with the small size of the system 
must introduce additional correlations which can modify 
the physic of the model. Indeed, in the present case, 
the maximum distance between two sites is about 
three elementary cells. The slices number is another 
problem. In principle, the method necessitates many 
simulations for increasing values of n with the aim of 
doing extrapolation. Simulations for only one value of n 
were realize because of the computation duration. Thus 
the transition temperatures should not be relevant. 
The choice of the sub-systems for the breakup of the 
system certainly has an influence. The breakup was 
chosen in such a way that the number of common sites 
of two neighbouring systems is minimum. This is in 
respect with the aim of the method. Moreover, this 
breakup is certainly the best in accordance with the 
physic of the Cu02 plane. The matrix size is not too 
large for the computation. 

Despite these weaknesses and the approximations, this 
simulation method is the alone method which allows 
to observe structures like loops and chains in the real 
two-dimensional array. It gives resuls which seem to 
relate the circulating current and the stripe concepts. 
Moreover the behaviour of the specific heat is in good 
qualitative agreement with the experimental data. Thus 
the presented simulation results can, perhaps, suggest 
a way to supplement the theories proposed to explain 
the hight Tc superconductivity. In this possible scenario 
the superconductivity is related to the TO chains, this 
means that the two holes of the copper end sites of each 
TO chain form a Cooper pair. The coherence length is 
the mean length between the two end sites of the TO 
chains. This mean length is about 2.7 times the lattice 
parameter (curve not shown) . This is in good agreement 
with the experimental data. The pairing process is the 
consequence of the particular structure of the Cu02 
plane and the repulsive interactions. Indeed, we observe 
(simulation results not shown) that the low temperature 
behaviour is strongly infuenced by the values of the 
interaction parameters Ud, Up and Vdp. For example, 
the low temperature metallic behaviour disappears when 
Vdp is larger than l.llSeV^, the other parameters being 
unchanged. This behaviour change is very sudden. 
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